Paso 3.2

Análisis de clases latentes: modelo seleccionado sin predictores, caracterización de clases y medidas de ajuste (glca)

Andrés González Santa Cruz
May 04, 2023
Show code
script src = "https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js"
Show code
 $(document).ready(function() {
    $('body').prepend('<div class=\"zoomDiv\"><img src=\"\" class=\"zoomImg\"></div>');
    // onClick function for all plots (img's)
    $('img:not(.zoomImg)').click(function() {
      $('.zoomImg').attr('src', $(this).attr('src')).css({width: '100%'});
      $('.zoomDiv').css({opacity: '1', width: 'auto', border: '1px solid white', borderRadius: '5px', position: 'fixed', top: '50%', left: '50%', marginRight: '-50%', transform: 'translate(-50%, -50%)', boxShadow: '0px 0px 50px #888888', zIndex: '50', overflow: 'auto', maxHeight: '100%'});
    });
    // onClick function for zoomImg
    $('img.zoomImg').click(function() {
      $('.zoomDiv').css({opacity: '0', width: '0%'}); 
    });
  });
  
Show code
<script src="hideOutput.js"></script> 
Show code
$(document).ready(function() {    
    $chunks = $('.fold');    
    $chunks.each(function () {      // add button to source code chunks     
    if ( $(this).hasClass('s') ) {       
        $('pre.r', this).prepend("<div class=\"showopt\">Show Source</div><br style=\"line-height:22px;\"/>");
            $('pre.r', this).children('code').attr('class', 'folded');     
            }      // add button to output chunks     
        if ( $(this).hasClass('o') ) {       
            $('pre:not(.r)', this).has('code').prepend("<div class=\"showopt\">Show Output</div><br style=\"line-height:22px;\"/>");       
            $('pre:not(.r)', this).children('code:not(r)').addClass('folded');        // add button to plots       
            $(this).find('img').wrap('<pre class=\"plot\"></pre>');       
            $('pre.plot', this).prepend("<div class=\"showopt\">Show Plot</div><br style=\"line-height:22px;\"/>");       
            $('pre.plot', this).children('img').addClass('folded');      
            }   
});    // hide all chunks when document is loaded   
    $('.folded').css('display', 'none')    // function to toggle the visibility   
    $('.showopt').click(function() {     
            var label = $(this).html();     
            if (label.indexOf("Show") >= 0) {       
                $(this).html(label.replace("Show", "Hide"));     
            } else {
              $(this).html(label.replace("Hide", "Show"));     
            }     
    $(this).siblings('code, img').slideToggle('fast', 'swing');   
    }); 
}); 

Cargamos los datos

Show code
rm(list = ls());gc()
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 535963 28.7    1209467 64.6   643711 34.4
Vcells 907206  7.0    8388608 64.0  1649632 12.6
Show code
#cargar glca sin resultado distal
load("data2_lca23_2023_04_29.RData")
#cargar glca con resultado distal
load("data2_lca2_adj4_2023_04_28.RData")
#data2_lca23.RData
#data2_lca2_adj4.RData

Cargamos los paquetes

Show code
knitr::opts_chunk$set(echo = TRUE)

if(!require(poLCA)){install.packages("poLCA")}
if(!require(poLCAParallel)){devtools::install_github("QMUL/poLCAParallel@package")}
if(!require(compareGroups)){install.packages("compareGroups")}
if(!require(parallel)){install.packages("parallel")}
if(!require(Hmisc)){install.packages("Hmisc")}
if(!require(tidyverse)){install.packages("tidyverse")}
try(if(!require(sjPlot)){install.packages("sjPlot")})
if(!require(emmeans)){install.packages("emmeans")}
if(!require(nnet)){install.packages("nnet")}
if(!require(here)){install.packages("here")}
if(!require(doParallel)){install.packages("doParallel")}
if(!require(progress)){install.packages("progress")}
if(!require(caret)){install.packages("caret")}
if(!require(rpart)){install.packages("rpart")}
if(!require(rpart.plot)){install.packages("rpart.plot")}
if(!require(partykit)){install.packages("partykit")}
if(!require(randomForest)){install.packages("randomForest")}
if(!require(ggcorrplot)){install.packages("ggcorrplot")}
if(!require(polycor)){install.packages("polycor")}
if(!require(tableone)){install.packages("tableone")}
if(!require(broom)){install.packages("broom")}
if(!require(plotly)){install.packages("plotly")}
if(!require(rsvg)){install.packages("rsvg")}
if(!require(DiagrammeRsvg)){install.packages("DiagrammeRsvg")}
if(!require(effects)){install.packages("effects")}
if(!require(glca)){install.packages("glca")}

Seleccionar modelos finales

Sin resultado distal

Show code
best_model_glca<- eval(parse(text = paste0("lca",best_model))) 
summary(best_model_glca)

Call:
glca(formula = f_preds, data = mydata_preds3, nclass = 7, n.init = 500, 
    decreasing = T, testiter = 500, maxiter = 10000, seed = seed, 
    verbose = FALSE)

Manifest items : CAUSAL EDAD_MUJER_REC PUEBLO_ORIGINARIO_REC PAIS_ORIGEN_REC HITO1_EDAD_GEST_SEM_REC MACROZONA PREV_TRAMO_REC 

Categories for manifest items :
                        Y = 1 Y = 2 Y = 3 Y = 4 Y = 5 Y = 6
CAUSAL                      2     3     4                  
EDAD_MUJER_REC              1     2     3     4     5     6
PUEBLO_ORIGINARIO_REC       1     2     3                  
PAIS_ORIGEN_REC             1     2     3                  
HITO1_EDAD_GEST_SEM_REC     1     2     3     4     5     6
MACROZONA                   1     2     3     4     5     6
PREV_TRAMO_REC              1     2     3     4     5      

Model : Latent class analysis 

Number of latent classes : 7 
Number of observations : 3789 
Number of parameters : 181 

log-likelihood : -26835.6 
     G-squared : 3742.786 
           AIC : 54033.2 
           BIC : 55162.61 

Marginal prevalences for latent classes :
Class 1 Class 2 Class 3 Class 4 Class 5 Class 6 Class 7 
0.19017 0.04980 0.09476 0.34128 0.12731 0.04374 0.15293 

Class prevalences by group :
    Class 1 Class 2 Class 3 Class 4 Class 5 Class 6 Class 7
ALL 0.19017  0.0498 0.09476 0.34128 0.12731 0.04374 0.15293

Item-response probabilities :
CAUSAL 
         Y = 1  Y = 2  Y = 3
Class 1 0.9688 0.0312 0.0000
Class 2 0.7077 0.2923 0.0000
Class 3 0.3454 0.6546 0.0000
Class 4 0.1273 0.8727 0.0000
Class 5 0.0756 0.9244 0.0000
Class 6 0.0556 0.0000 0.9444
Class 7 0.0086 0.0000 0.9914
EDAD_MUJER_REC 
         Y = 1  Y = 2  Y = 3  Y = 4  Y = 5  Y = 6
Class 1 0.0006 0.0106 0.2417 0.4167 0.2945 0.0359
Class 2 0.0684 0.0360 0.2950 0.3044 0.2476 0.0486
Class 3 0.0000 0.0170 0.2900 0.3167 0.3140 0.0623
Class 4 0.0017 0.0414 0.2874 0.3261 0.2783 0.0651
Class 5 0.0030 0.0017 0.0870 0.2979 0.5370 0.0733
Class 6 0.0000 0.1794 0.3665 0.2473 0.1770 0.0299
Class 7 0.0017 0.3681 0.2947 0.2092 0.1071 0.0192
PUEBLO_ORIGINARIO_REC 
         Y = 1  Y = 2  Y = 3
Class 1 0.0859 0.8644 0.0496
Class 2 1.0000 0.0000 0.0000
Class 3 0.1255 0.8373 0.0372
Class 4 0.1239 0.8184 0.0577
Class 5 0.1401 0.8599 0.0000
Class 6 0.1290 0.8337 0.0373
Class 7 0.1521 0.8030 0.0449
PAIS_ORIGEN_REC 
         Y = 1  Y = 2  Y = 3
Class 1 0.0000 0.8801 0.1199
Class 2 0.0838 0.8285 0.0876
Class 3 0.0061 0.0000 0.9939
Class 4 0.0000 0.9887 0.0113
Class 5 0.0000 0.9268 0.0732
Class 6 0.0000 0.0000 1.0000
Class 7 0.0000 0.9922 0.0078
HITO1_EDAD_GEST_SEM_REC 
         Y = 1  Y = 2  Y = 3  Y = 4  Y = 5  Y = 6
Class 1 0.0754 0.2027 0.2229 0.4909 0.0000 0.0081
Class 2 0.0500 0.0000 0.0914 0.2432 0.3096 0.3058
Class 3 0.0077 0.0071 0.3264 0.3672 0.2250 0.0666
Class 4 0.0063 0.0000 0.3262 0.3622 0.2232 0.0821
Class 5 0.0080 0.0040 0.6071 0.2574 0.1040 0.0195
Class 6 0.0202 0.7950 0.1848 0.0000 0.0000 0.0000
Class 7 0.0088 0.7795 0.2117 0.0000 0.0000 0.0000
MACROZONA 
         Y = 1  Y = 2  Y = 3  Y = 4  Y = 5  Y = 6
Class 1 0.0000 0.3995 0.1653 0.1662 0.1128 0.1562
Class 2 0.0311 0.1991 0.3146 0.1668 0.0996 0.1888
Class 3 0.0000 0.6108 0.0844 0.0372 0.2360 0.0316
Class 4 0.0009 0.3192 0.1827 0.2115 0.0888 0.1970
Class 5 0.0000 0.6956 0.1042 0.0629 0.0663 0.0711
Class 6 0.0060 0.5570 0.0803 0.0159 0.2994 0.0413
Class 7 0.0052 0.3843 0.1750 0.1448 0.0859 0.2048
PREV_TRAMO_REC 
         Y = 1  Y = 2  Y = 3  Y = 4  Y = 5
Class 1 0.0000 0.0915 0.5565 0.3520 0.0000
Class 2 0.0276 0.0370 0.6457 0.2347 0.0550
Class 3 0.0048 0.0000 0.6095 0.2880 0.0977
Class 4 0.0008 0.0032 0.6439 0.3521 0.0000
Class 5 0.0000 0.7662 0.0349 0.1908 0.0082
Class 6 0.0300 0.0043 0.5214 0.1636 0.2808
Class 7 0.0017 0.0701 0.7223 0.2006 0.0053
Show code
#https://rdrr.io/cran/glca/src/R/summary.glca.R
#Class prevalence by group mean(best_model_glca_w_distal_outcome$posterior$ALL$`Class 1`)
#Item-response probabilities (most likely response)  print(sapply(best_model_glca_w_distal_outcome$param$rho[[1]],
#print(sapply(best_model_glca_w_distal_outcome$param$rho$ALL, function(m) apply(m, 1, which.max)))
#function(m) apply(m, 1, which.max)))
Show code
#https://rdrr.io/cran/glca/src/R/plot.glca.R
plot(eval(parse(text = paste0("lca",best_model))), ask=F)
Selected Model

Figure 1: Selected Model

Selected Model

Figure 2: Selected Model

Selected Model

Figure 3: Selected Model

Selected Model

Figure 4: Selected Model

Selected Model

Figure 5: Selected Model

Selected Model

Figure 6: Selected Model

Selected Model

Figure 7: Selected Model

Selected Model

Figure 8: Selected Model

Show code
rho_glca<- 
do.call("bind_rows",best_model_glca$param$rho$ALL) %>% 
  t() %>% 
  round(2) %>% 
  data.table::data.table(keep.rownames = T) %>% 
  magrittr::set_colnames(c("variables", paste0("Class",1:7))) %>% 
  tidyr::separate(variables, into=c("var", "prob"), sep=".Y =")

lcmodel_glca <- reshape2::melt(rho_glca, level=2) %>% dplyr::rename("class"="variable")

traductor_cats <- readxl::read_excel("tabla12.xlsx") %>% 
  dplyr::mutate(level=readr::parse_double(level)) %>% 
  dplyr::mutate(level= dplyr::case_when(grepl("CAUSAL",charactersitic)~ level-1,T~level)) %>% 
  dplyr::mutate(charactersitic=gsub(" \\(%\\)", "", charactersitic))



lcmodel_glca<- lcmodel_glca %>% 
  dplyr::mutate(pr=as.numeric(gsub("[^0-9.]+", "", prob))) %>% 
  dplyr::left_join(traductor_cats[,c("charactersitic", "level", "CATEGORIA")], by= c("var"="charactersitic", "pr"="level"))

lcmodel_glca$text_label<-paste0("Categoria:",lcmodel_glca$CATEGORIA,"<br>%: ",scales::percent(lcmodel_glca$value))

zp3 <- ggplot(lcmodel_glca,aes(x = var, y = value, fill = factor(pr), label=text_label))
zp3 <- zp3 + geom_bar(stat = "identity", position = "stack")
zp3 <- zp3 + facet_grid(class ~ .) 
zp3 <- zp3 + scale_fill_brewer(type="seq", palette="Greys", na.value = "white") +theme_bw()
zp3 <- zp3 + labs(y = "Porcentaje de probabilidad de respuesta", 
                  x = "",
                  fill ="Cateorías de\nRespuesta")
zp3 <- zp3 + theme( axis.text.y=element_blank(),
                    axis.ticks.y=element_blank(),                    
                    panel.grid.major.y=element_blank())
zp3 <- zp3 + guides(fill = guide_legend(reverse=TRUE))
zp3 <- zp3 + theme(axis.text.x = element_text(angle = 30, hjust = 1))
#print(zp1)

ggplotly(zp3, tooltip = c("text_label"))%>% layout(xaxis= list(showticklabels = T),height=600, width=800)

Figure 9: Selected Model

Show code
ggsave("_fig3_LCA_distribuciones_glca.png",zp3, dpi= 600)

lcmodel_glca %>% rio::export("variables_probabilities_in_category_glca.xlsx")
Error in saveWorkbook(wb, file = file, overwrite = overwrite): File already exists!


Con resultado distal

Show code
best_model_glca_w_distal_outcome<- eval(parse(text = paste0("lca2",best_model2)))

summary(best_model_glca_w_distal_outcome)

Call:
glca(formula = f_adj, data = mydata_preds3, nclass = 7, n.init = 500, 
    decreasing = T, testiter = 500, maxiter = 10000, seed = seed, 
    verbose = FALSE)

Manifest items : CAUSAL EDAD_MUJER_REC PUEBLO_ORIGINARIO_REC PAIS_ORIGEN_REC HITO1_EDAD_GEST_SEM_REC MACROZONA PREV_TRAMO_REC 
Covariates (Level 1) : outcome 

Categories for manifest items :
                        Y = 1 Y = 2 Y = 3 Y = 4 Y = 5 Y = 6
CAUSAL                      2     3     4                  
EDAD_MUJER_REC              1     2     3     4     5     6
PUEBLO_ORIGINARIO_REC       1     2     3                  
PAIS_ORIGEN_REC             1     2     3                  
HITO1_EDAD_GEST_SEM_REC     1     2     3     4     5     6
MACROZONA                   1     2     3     4     5     6
PREV_TRAMO_REC              1     2     3     4     5      

Model : Latent class analysis 

Number of latent classes : 7 
Number of observations : 3789 
Number of parameters : 187 

log-likelihood : -26763.76 
     G-squared : 5115.172 
           AIC : 53901.52 
           BIC : 55068.38 

Marginal prevalences for latent classes :
Class 1 Class 2 Class 3 Class 4 Class 5 Class 6 Class 7 
0.18988 0.04756 0.08593 0.31208 0.16780 0.04381 0.15293 

Class prevalences by group :
    Class 1 Class 2 Class 3 Class 4 Class 5 Class 6 Class 7
ALL 0.18988 0.04756 0.08593 0.31208  0.1678 0.04381 0.15293

Logistic regression coefficients :
            Class 1/7 Class 2/7 Class 3/7 Class 4/7 Class 5/7
(Intercept)    0.9933   -1.0814   -0.0580    1.9176   -0.3817
outcome1      -0.8885   -0.0947   -0.5814   -1.4400    0.5080
            Class 6/7
(Intercept)   -1.5881
outcome1       0.3634
Item-response probabilities :
CAUSAL 
         Y = 1  Y = 2  Y = 3
Class 1 0.9716 0.0284 0.0000
Class 2 0.7457 0.2543 0.0000
Class 3 0.3733 0.6267 0.0000
Class 4 0.1303 0.8697 0.0000
Class 5 0.0748 0.9252 0.0000
Class 6 0.0591 0.0000 0.9409
Class 7 0.0080 0.0000 0.9920
EDAD_MUJER_REC 
         Y = 1  Y = 2  Y = 3  Y = 4  Y = 5  Y = 6
Class 1 0.0007 0.0112 0.2433 0.4167 0.2923 0.0359
Class 2 0.0681 0.0362 0.2959 0.3122 0.2496 0.0380
Class 3 0.0000 0.0169 0.3047 0.3113 0.3103 0.0569
Class 4 0.0017 0.0464 0.3004 0.3202 0.2646 0.0666
Class 5 0.0034 0.0000 0.1020 0.3156 0.5052 0.0738
Class 6 0.0000 0.1787 0.3661 0.2472 0.1782 0.0298
Class 7 0.0017 0.3683 0.2944 0.2092 0.1071 0.0192
PUEBLO_ORIGINARIO_REC 
         Y = 1  Y = 2  Y = 3
Class 1 0.0844 0.8658 0.0498
Class 2 1.0000 0.0000 0.0000
Class 3 0.1309 0.8321 0.0369
Class 4 0.1286 0.8071 0.0642
Class 5 0.1382 0.8618 0.0000
Class 6 0.1292 0.8336 0.0372
Class 7 0.1519 0.8031 0.0449
PAIS_ORIGEN_REC 
         Y = 1  Y = 2  Y = 3
Class 1 0.0000 0.8834 0.1166
Class 2 0.0848 0.8242 0.0910
Class 3 0.0084 0.0000 0.9916
Class 4 0.0000 0.9692 0.0308
Class 5 0.0000 0.9219 0.0781
Class 6 0.0000 0.0000 1.0000
Class 7 0.0000 0.9917 0.0083
HITO1_EDAD_GEST_SEM_REC 
         Y = 1  Y = 2  Y = 3  Y = 4  Y = 5  Y = 6
Class 1 0.0749 0.2015 0.2306 0.4870 0.0000 0.0059
Class 2 0.0574 0.0000 0.0905 0.2468 0.2951 0.3101
Class 3 0.0077 0.0077 0.3147 0.3899 0.2207 0.0594
Class 4 0.0013 0.0000 0.2861 0.3710 0.2451 0.0964
Class 5 0.0161 0.0047 0.6078 0.2572 0.0992 0.0149
Class 6 0.0204 0.7946 0.1851 0.0000 0.0000 0.0000
Class 7 0.0088 0.7794 0.2119 0.0000 0.0000 0.0000
MACROZONA 
         Y = 1  Y = 2  Y = 3  Y = 4  Y = 5  Y = 6
Class 1 0.0000 0.3971 0.1643 0.1693 0.1109 0.1584
Class 2 0.0327 0.1899 0.3338 0.1624 0.0986 0.1826
Class 3 0.0000 0.6294 0.0846 0.0149 0.2518 0.0193
Class 4 0.0009 0.2832 0.1819 0.2321 0.0903 0.2116
Class 5 0.0000 0.6807 0.1169 0.0600 0.0713 0.0712
Class 6 0.0060 0.5588 0.0796 0.0161 0.2983 0.0412
Class 7 0.0052 0.3843 0.1752 0.1447 0.0860 0.2046
PREV_TRAMO_REC 
         Y = 1  Y = 2  Y = 3  Y = 4  Y = 5
Class 1 0.0000 0.0874 0.5601 0.3524 0.0000
Class 2 0.0243 0.0403 0.6488 0.2295 0.0571
Class 3 0.0066 0.0000 0.6076 0.2792 0.1065
Class 4 0.0013 0.0116 0.6462 0.3410 0.0000
Class 5 0.0000 0.5706 0.1718 0.2509 0.0068
Class 6 0.0299 0.0040 0.5209 0.1640 0.2812
Class 7 0.0017 0.0703 0.7221 0.2006 0.0053
Show code
rho_glca_adj<- 
do.call("bind_rows",best_model_glca_w_distal_outcome$param$rho$ALL) %>% 
  t() %>% 
  round(2) %>% 
  data.table::data.table(keep.rownames = T) %>% 
  magrittr::set_colnames(c("variables", paste0("Class",1:7))) %>% 
  tidyr::separate(variables, into=c("var", "prob"), sep=".Y =")

lcmodel_glca_adj <- reshape2::melt(rho_glca_adj, level=2) %>% dplyr::rename("class"="variable")

lcmodel_glca_adj<- lcmodel_glca_adj %>% 
  dplyr::mutate(pr=as.numeric(gsub("[^0-9.]+", "", prob))) %>% 
  dplyr::left_join(traductor_cats[,c("charactersitic", "level", "CATEGORIA")], by= c("var"="charactersitic", "pr"="level"))

lcmodel_glca_adj$text_label<-paste0("Categoria:",lcmodel_glca_adj$CATEGORIA,"<br>%: ",scales::percent(lcmodel_glca_adj$value))

zp4 <- ggplot(lcmodel_glca_adj,aes(x = var, y = value, fill = factor(pr), label=text_label))
zp4 <- zp4 + geom_bar(stat = "identity", position = "stack")
zp4 <- zp4 + facet_grid(class ~ .) 
zp4 <- zp4 + scale_fill_brewer(type="seq", palette="Greys", na.value = "white") +theme_bw()
zp4 <- zp4 + labs(y = "Porcentaje de probabilidad de respuesta", 
                  x = "",
                  fill ="Cateorías de\nRespuesta")
zp4 <- zp4 + theme( axis.text.y=element_blank(),
                    axis.ticks.y=element_blank(),                    
                    panel.grid.major.y=element_blank())
zp4 <- zp4 + guides(fill = guide_legend(reverse=TRUE))
zp4 <- zp4 + theme(axis.text.x = element_text(angle = 30, hjust = 1))
#print(zp1)

ggplotly(zp4, tooltip = c("text_label"))%>% layout(xaxis= list(showticklabels = T),height=600, width=800)

Figure 10: Selected Model

Show code
ggsave("_fig4_LCA_distribuciones_glca_adj.png",zp4, dpi= 600)

lcmodel_glca %>% rio::export("variables_probabilities_in_category_glca_adj.xlsx")
Error in saveWorkbook(wb, file = file, overwrite = overwrite): File already exists!
Show code
#https://rdrr.io/cran/glca/src/R/plot.glca.R
plot(eval(parse(text = paste0("lca2",best_model2))), ask=F)
Selected Model

Figure 11: Selected Model

Selected Model

Figure 12: Selected Model

Selected Model

Figure 13: Selected Model

Selected Model

Figure 14: Selected Model

Selected Model

Figure 15: Selected Model

Selected Model

Figure 16: Selected Model

Selected Model

Figure 17: Selected Model

Selected Model

Figure 18: Selected Model

Show code
glca_gam_dist_outcome<- best_model_glca_w_distal_outcome$param$gamma[[1]]
glca_gam_dist_outcome<-colnames(paste0("Class",1:7))
#conditional probabilities
#Pr(B1=1|Class 3)
posteriors_glca_adj <- data.frame(best_model_glca_w_distal_outcome$posterior, 
                             predclass=best_model_glca_w_distal_outcome$param$gamma) 

classification_table_adj <- plyr::ddply(posteriors, "predclass", function(x) colSums(x[,1:length(LCA_best_model_adj_mod$P)]))
clasification_errors_adj<- 1 - sum(diag(as.matrix(classification_table_adj[,2:(length(LCA_best_model_adj_mod$P)+1)]))) / sum(classification_table_adj[,2:(length(LCA_best_model_adj_mod$P)+1)]) 

warning(paste("Error de clasificación: ", round(clasification_errors_adj,2)))


entropy_alt <- function(p) sum(-p * log(p))
error_prior_adj <- entropy_alt(LCA_best_model_adj_mod$P) # Class proportions
error_post_adj <- mean(apply(LCA_best_model_adj_mod$posterior, 1, entropy_alt),na.rm=T)
R2_entropy_alt_adj <- (error_prior_adj - error_post_adj) / error_prior_adj
warning(paste("Entropía: ", round(R2_entropy_alt_adj,2)))


#https://stackoverflow.com/questions/72783185/entropy-calculation-gives-nan-is-applying-na-omit-a-valid-tweak
entropy.safe <- function (p) {
  if (any(p > 1 | p < 0)) stop("probability must be between 0 and 1")
  log.p <- numeric(length(p))
  safe <- p != 0
  log.p[safe] <- log(p[safe])
  sum(-p * log.p)
}

error_prior2_adj <- entropy.safe(LCA_best_model_adj_mod$P) # Class proportions
error_post2_adj <- mean(apply(LCA_best_model_adj_mod$posterior, 1, entropy.safe),na.rm=T)
R2_entropy_safe_adj <- (error_prior2_adj - error_post2_adj) / error_prior2_adj
warning(paste("Entropía segura: ", round(R2_entropy_safe,2)))

#https://gist.github.com/daob/c2b6d83815ddd57cde3cebfdc2c267b3
warning(paste("Entropía (solución Oberski): ", round(entropy.R2(LCA_best_model_adj_mod),2)))

#\#minimum average posterior robability of cluster membership (\>0.7), interpretability (classes are clearly distinguishable), and parsimony (each class has a sufficient sample size for further analysis; n≥50).


Show code
coef(eval(parse(text = paste0("lca2",best_model2))))
Class 1 / 7 :
            Odds Ratio Coefficient  Std. Error  t value  Pr(>|t|)   
(Intercept)     0.2993     -1.2064      0.4550   -2.651   0.00805 **
outcome1        0.8654     -0.1445      0.4527   -0.319   0.74957   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Class 2 / 7 :
            Odds Ratio Coefficient  Std. Error  t value  Pr(>|t|)  
(Intercept)     1.4648      0.3817      0.3159    1.208    0.2271  
outcome1        0.6017     -0.5080      0.3070   -1.654    0.0981 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Class 3 / 7 :
            Odds Ratio Coefficient  Std. Error  t value  Pr(>|t|)    
(Intercept)     9.9678      2.2994      0.2976    7.726  1.43e-14 ***
outcome1        0.1426     -1.9480      0.2784   -6.997  3.11e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Class 4 / 7 :
            Odds Ratio Coefficient  Std. Error  t value  Pr(>|t|)    
(Intercept)     1.3822      0.3237      0.3406    0.950  0.341966    
outcome1        0.3364     -1.0894      0.3303   -3.298  0.000983 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Class 5 / 7 :
            Odds Ratio Coefficient  Std. Error  t value  Pr(>|t|)    
(Intercept)     3.9550      1.3750      0.2990    4.599  4.39e-06 ***
outcome1        0.2475     -1.3964      0.2854   -4.893  1.04e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Class 6 / 7 :
            Odds Ratio Coefficient  Std. Error  t value  Pr(>|t|)
(Intercept)     0.4967     -0.6997      0.5021   -1.394     0.164
outcome1        0.5473     -0.6027      0.4972   -1.212     0.226
Show code
save.image("data2_lca3_glca.RData")
Show code
require(tidyverse)
sesion_info <- devtools::session_info()
dplyr::select(
  tibble::as_tibble(sesion_info$packages),
  c(package, loadedversion, source)
) %>% 
  DT::datatable(filter = 'top', colnames = c('Row number' =1,'Variable' = 2, 'Percentage'= 3),
              caption = htmltools::tags$caption(
        style = 'caption-side: top; text-align: left;',
        '', htmltools::em('Packages')),
      options=list(
initComplete = htmlwidgets::JS(
        "function(settings, json) {",
        "$(this.api().tables().body()).css({
            'font-family': 'Helvetica Neue',
            'font-size': '50%', 
            'code-inline-font-size': '15%', 
            'white-space': 'nowrap',
            'line-height': '0.75em',
            'min-height': '0.5em'
            });",#;
        "}")))